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"^ ■ This paper describes the design and work of a set of computer routines capable for nu- 

C^^ , merical computation of generahzed parton distributions (GPDs) evolution at the leading 

order. The main intention of this work is to present a fast-working computer code making 
possible fitting of GPDs parameters to the data on hard electron-nucleon scattering. 
> 

(M ■ 1 Introduction 

^ I Generalized parton distributions EJ El IH E] parameterize matrix elements of quark and gluon 

O ' twist-2 operators between nucleon states with non-equal momenta. They can be represented as 

r-| ■ functions of fractions x and C, of longitudinal momentum of the proton carried by the parton 

^. and of t = (pi —p2)^, the square of the 4-momentum transfer between initial and final protons^ 

Ph! (see. Fig.P). 

> 
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Figure 1: Partonic picture for GPDs. The value of ^ is given by ^ = Xb/(2 — x^), where Xb is 
Bjorken variable. 

GPDs are the basis for description of hard elastic lepton-nucleon scattering in terms of 
QCD, just like the ordinary PDFs are the basis of QCD description of the DIS. Besides the 
description of hard elastic reactions, GPDs give access to the orbital momentum of partons in 
the nucleon [H] and 3-dimensional nucleon structure [7111113 El- The importance of GPDs for 
understanding of the nucleon structure raised a problem of their extraction from hard elastic 
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^The variable t does not enter the evolution equations and is therefore omitted in this paper. 



scattering data, first of all from deeply virtual Compton scattering (DVCS). The activity of 
measuring DVCS observables at lepton-nucleon facilities is quite intense [TTJ [T21 HSl UH Hi] • 

Similar to PDFs, GPDs are subject of perturbative QCD evolution which should be taken 
care of in extraction of GPDs from the data on lepton-nucleon scattering. The LO kernels for 
GPDs evolution are well-known [U El 01 IH El CEl IIZ| • The NLO results have been pubhshed jIHI. 
At present, there exists a code [iniEIII which can be used for leading and next-to-leading order 
evolution of GPDs. However, on a 1.7 GHz Pentium-4 machine this code takes about 15 seconds 
for LO evolution and about 40 seconds for the NLO which is not fast enough for performing 
fitting procedure requiring very large number of evolution routine calls. Evolution equations for 
GPDs at the LO can be also solved by a decomposition of GPDs in series of special functions 
PH 122] • The advantage of this method is that it does not require introduction of a difference 
scheme for Q^. However, the decomposition is not trivial for numerical implementation and the 
number of terms in the series required to achieve a reliable precision can be large (depending 
on the input GPDs). 

In the present paper another code is introduced^, which does the LO GPD evolution with 
a reasonable precision for ~ 10~^ seconds, making possible its usage for fitting. The fast work 
was achieved by optimization of the computational procedure. It uses 4-th order Runge-Kutta 
method for solving the evolution equations. Since in the broad DGLAP region (|x| > ^) the 
variation of GPDs is not sizable, and in the narrow ERBL region (|a:| < ^) is quite significant, 
the convolution integrals are performed on a logarithmic grid so that the number of points 
in ERBL and DGLAP regions is the same. The code itself is also optimized. For example, 
the operation of division is avoided as much as possible. This was done since division on the 
machines commonly used (such as Pentium- IV) is performed 6 times slower than multiplication. 
So, if an expression to be computed reads a/b/c, its optimized version is a/{h * c). Since the 
NLO expressions are considerably larger then the LO ones, their implementation is not done 
at the moment. 

2 Evolution kernels 

The LO kernels for GPDs evolution are well-known. In the code those published in [3 QEl ^| 
are used^. The numerical computation of the evolution requires that singularities in the kernel 
to be cancelled analytically. The kernels with regularized singularities in the form they take in 
the code are given in the Appendix. 

The kernels assume the definition of the GPDs given in J2S] • In the forward limit they obey 



Hns{x, ^ = 0, t = 0) = u{x) — d{x) , 
Hq{x, ^ = 0, t = 0) = u{x) + d{x) + s{x) , 

HG{x,^ = 0,t = 0)=xg{x), (1) 

Hq{x, ^ = 0, t = 0) = Au{x) + Ad{x) + As{x) , 
HG{x,^ = 0,t = 0)=xAg{x), 

withg(x) = qv{x)+q{x), Aq{x) = Aqy{x)+Aq{x), q{—x) = — g(x), g{—x) = —g{x), Aq{—x) = 
Aq{x), Ag{—x) = Ag{x). Note that quark singlet and non-singlet densities are not divided by 
Uf. This is important in the singlet sector, where non-diagonal parts of the kernel depend on the 



^The code is available at http://www.i&.de/~vinnikov/gpd_evol.tar.gz 

^Literally, the expressions from Ref. 17 were taken. Due to the symmetry relations of singlet GPDs, these 



kernels are equivalent to those published in Ref. [S], (see Eq. 5.21 of Ref. |S]). 



normalization of the quark singlet distribution. For non-singlet channel, isovector combination 
was chosen to represent the work of the code. It is clear that any other non-singlet combination 
of GPDs can be taken as well. 

In the singlet sector, the symmetry properties Hq^—x.^^t) = —HQ{x,^,t), Hc{—x,^,t) = 
Hc{—x, ^, t), Hq{—x, ^, t) = Hq{x, ^, t) and Hg{—x, C,, t) = —Hg{x, ^, t) are used. This allowed 
to reduce the computational time by a factor of 4. Thus, the singlet kernels are given only for 
x>0. 

For all parts of the kernel, convolution at a; = ±i^ was done separately. The kernels at 
X = ±^ were taken as x — > ^ limits of the corresponding DGLAP or (equivalently for x — > ±^) 
ERBL expressions. 

The running strong coupling constant is taken as a solution of the NLO RG equation 

1^ = -boalif^) - b,alif^) - hatif^) (2) 

where the coefficients bo and 6i are given by 

^ (33 -2A^/ 



f)^ 



VItx 
bi = ^(153-19iV^), (3) 

^2 = -TTTT^ (77139 -l5099Nf + 325N^f 

Here Np is the number of active quark flavors. The initial condition for Eq. |2l is taken as 
the value of as{Mz)- Then the equation is solved back to smaller values of fi. Down to 
fi = m{b) = 4.3 GeV the number of active flavors is taken to be equal 5. Between /x = m{c) = 1.3 
GeV and fi = in{b), Nf = 4 is used. Below /i = 'm(c) Nf is 3. Note that the number of active 
flavors Nf entering the ^^(/i) RG equation is different from the number of flavors in the proton 
Uf which enters into the quark-gluon evolution kernel. In the code, Uf is always equal to 3. 

3 Computation of the convolution integrals 

The integrals are computed on a logarithmic grid. The required parameterization for non-singlet 
case reads: 

Xi = -5e-(^-2") + 5 (-l<x<0, 0<i<2n), (4) 

Xi = 5e(*-2") -5 {0<x<l,2n<i<4:n), (5) 

where 7 and 5 are given by 

■> = lH\ + y»e-12i^4), *=^. (6) 

This grid has the following useful properties: 

^0 = ~1; Xn = — ^, X2n = 0, X^n = S^j X^n = 1 ( ' ) 

i.e. it provides the same level of detalization for ERBL and DGLAP regions. 

For the singlet case the grid size is reduced so that to cover only < x < 1 interval: 

Xi = 5e'' -5 (0 < i < 2n) , (8) 



so that 

xo = 0, x„ = ^, a;2„ = 1 . (9) 

The parameters 7 and 6 are the same as in the non-singlet case (Eq. ^. 

After introduction of the logarithmic grid the convolution integrals can be computed using 
the standard equidistant Simpson's formula taking into account the Jacobian for the grid: 

1 2n 4n 

f{x)dx = I i{5 - Xi)f{xi)di + hy{6 + Xi)f{xi)di . (10) 

-1 2n 

The ^ uncertainty in the integrals can be resolved with the L'Hopital's rule. Since for the 
numerical integration the Simpson's {0{Ax^)) method is used, the numerical computation 
of the derivative requires 0{Ax^) precision. Depending on the position of the singularity 
(availability of other grid points on the left and on the right), one of the following formulas is 
be used: 

^'(^*) = TTTT—^ (-ll/(^^) + 18/(a:m) - 9/(a;,+2) + 2f{x,,+,)) , (11) 

07(0 + X) 

/'(^^) = « ^L , (-2/(a;.-i) - Sfix,) + 6/(a;,+i) - /(x,+2)) , (12) 

67(d) + X) 

/'(^^) = TTTT—-^ (^(^-2) - 6/(x,_i) + 3/(a;,) + 2f{x,+,)) , (13) 

07(0 + X) 

67(0 + X) 

For a; < the denominator should be replaced by 67(5 — x). 

To verify the convolution code, a computation was done for a simple functional form GPDs. 
For the non-singlet, spin- independent quark-gluon and gluon-gluon parts and spin-dependent 
quark-quark and gluon-quark parts 

Ht,,tix,^,Q^) = l-x^ (15) 

was taken, and for singlet spin-independent quark-quark and gluon-quark parts and spin- 
dependent quark-gluon and gluon-gluon parts 

Htestix,^,Q^)=x-x'. (16) 

It was found that taking n = 50 provides 4-6 digits precision for all the integrals. It should 
be noted that the test functions are not convenient for the Simpson method calculation since 
in the logarithmic grid they fall steeply when x — >■ 1. However, they made possible analytical 
computation of the integrals and therefore they were used in the test. The real GPDs which 
fall rapidly as x — ^ ±1 should be computed with use of Simpson's method more precisely on 
such a grid. 

4 The evolution procedure 

Basing on the integrals obtained in the previous section, the code for GPD evolution was 
written. It is based on the 4-th order Runge-Kutta method. Indeed, after introduction of the 
grid (0) the convolution integrals are transformed into sums. Therefore, the integro-differential 

4 



evolution equation is transformed into a set of 4n + 1 ordinary differential equations of the 
form^ 

^^ = F,(lnQ^ Ho, H,, ...H,,,), (17) 

where Hi = H{xi,C,,t,\iaQ). The 4-th order Runge-Kutta scheme on the InQ^-grid with m 
nodes 

therefore reads: 
where^ 



IngJ = \nQl + jd\nQ^, < j < m - 1 (18) 

Hi = Hr + ^ (f^^ + 2kP + 2kP + kf) , (19) 



ki'^=FSnQlHo,H,,...H,^), (20) 



kr - F.dnQ? + ^, f^o + '-^ki^\H, + ^.«, ...H,. + ^fcS, ), (21) 

k?^ - FSn Q? + ^, i^o + '-^kl^\H, + ^M^ ...i^.. + ^€\ ), (22) 
fcf ) = Fi(lnQ2 + d InQ^, Hq + d InQ^ k^^) ^ H^ + d InQ' fcf \ ...i74n + d\nQ^ k^!,) • (23) 

The results of the computation are presented in Figs.EHZIfor typical kinematics of HERMES 
and HERA. The code works both in forward (from smaller to larger Q"^) and backward (from 
larger to smaller Q"^) directions. As it can be seen, choice of ra = 20 and m = 1 already gives 
satisfactory results. Distinctions from the careful solution with n = 500 and m = 500 can 
be noticed around x = ^, but they are much smaller then typical experimental errors, so the 
choice of n = 20, m = 1 can be recommended for time-consuming fitting procedures. If the 
time limitations are not very strict, larger number of n and m can be taken providing higher 
precision. 

The evolution procedure was checked for satisfying the momentum conservation sum rule 



d 



dlnQ2 

-1 

and flavor conservation sum rule 

d 



1 

J {xHQix,tQ') + HGix,tQ'))dx = (24) 



dlnQ 



1 

J HNs{x,^,Q^)dx = 0. (25) 



-1 



For n = 20, m = 1, the sum rules are satisfied within 0.01% accuracy at HERMES kinematics 
and 0.5% accuracy at HERA kinematics. 

Computational time depends on n and m as 

t = an^m, (26) 

where m is the number of intervals for In Q. The coefficient a depends on the computational 
power of the machine. For Pentium-IV 1.7 GHz a ~ 1.4 ■ 10~^ sec in the non-singlet case and ~ 
2 ■ 10~^ in the singlet case. This means that computation of evolution with n=20, m=l takes 
about 10"'^ sec. 



*In singlet case, the function index runs from to 2n instead of to An in the non-singlet sector. 
^R.h.s. of these equations depend explicitly (i.e. not via GPDs) on Q^ through as{Q^) 
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Figure 2: Evolution of the quark non-singlet GPD H^s for ^ = 0.1 from Q^ = 4 GeV^ to 1 
GeV2 (left panel) and for ^ = 10"^ from Q^ = A GeV^ to 40 GeV^ (right panel). The initial 
GPD is obtained using the double-distribution ansatz j21l 123 with 6 = 1 and CTEQ6M 120] 
PDF at the scale of 4 GeV^. For the evolved distribution, n and m represent parameters of the 
numerical evolution procedure described in Eqs. (0) and p8|) correspondingly. 
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Figure 3: Evolution of the quark singlet GPD Hq for ^ = 0.1 from Q^ = 4 GeV^ to 1 GeV^ 
(left panel) and for ^ = 10~^ from Q^ = 4 GeV^ to 40 GeV^ (right panel). See caption of Fig.|21 
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Figure 4: Evolution of the gluon GPD Hq for ^ = 0.1 from Q^ = A GeV^ to 1 GeV^ (left panel) 
and for ^ = lO^'* from Q'^ = A GeV^ to 40 GeV^ (right panel). See caption of Fig. |21 (the input 
profile parameter 6 = 2 was taken). 
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Figure 5: Evolution of the spin-dependent non-singlet GPD Hjsrs for 'C = 0-1 from Q^ = 4 GeV^ 
to 1 GeV2 (left panel) and for ^ = 10"^ from Q^ = 4 QeV^ to 40 GeV^ (right panel). The 
initial GPD is obtained using the double-distribution ansatz [2^ |23] with 6=1 and Bliimlein- 
Bottcher j27j PDF at the scale of 4 GeV^. For the evolved distribution, n and m represent 
parameters of the numerical evolution procedure described in Eqs. and ^1 correspondingly. 
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Evolution of the spin-dependent singlet quark GPD Hq for ,^ = 0.1 from Q^ = 4 
GeV^ (left panel) and for ^ = 10"^ from Q^ = 4 GeV^ to 40 GeV^ (right panel). See 
Fig. El 
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Figure 7: Evolution of the spin-dependent gluon GPD Hq for ^ = 0.1 from Q"^ = A GeV^ to 1 
GeV^ (left panel) and for ^ = lO"'* from Q'^ = A GeV^ to 40 GeV^ (right panel). See caption 
of Fig. El (the input profile parameter 6 = 2 was taken). 
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Figure 8: Influence of a small (0.003 of the corresponding GPD value) hand-made spike at 
X = 0.105 in the input GPD with ^ = 0.1 on the output quark (left) and gluon (right) GPDs. 

5 Existing numerical instabilities 

Even though the code is quite effective, instabilities stemming from the complicate structure 
of the kernels still arise in the output. So far, only one instability was found. It is related to 
the fact that due to the "+" regularization procedure the convolution integral at a given value 
of X depends to a large extent on the value of derivative of the GPD at this point. Therefore, 
if an input GPD has sudden spikes (not related to a possible semi-singular behavior at x = ^ 
which the code can handle), the code can amplify them. This type of instability is illustrated in 
Fig. IHl The problem occurs only in the gluon-gluon part of the kernel and only for evolution in 
the backward direction, i.e. from larger Q^ to smaller one. It should be noted that the smaller 
is the step in x, the bigger is the value of the derivative the sudden spikes produce. Therefore, 
when for some reason the task is to perform evolution with very high precision, the stability of 
the input GPDs must be taken care of. 

6 Summary 

A code for prompt numerical computation of leading order GPDs evolution is presented. Its 
fundamentals are computation of the convolution integrals on a logarithmic grid and 4-th order 
Runge-Kutta method for integration of the differential equations. The code makes possible 
fast computation of the GPD evolution. On a 1.7 GHz Pentium-4 machine computational time 
required to perform evolution in the Q^ range from 4 to 1 GeV^ or from 4 to 40 GeV^ is of 
order of 10~^ seconds. Such code makes possible fitting of GPDs parameters from data on hard 
elastic lepton-nucleon reactions. 
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A Explicit expressions for the kernels 

A.l Non-singlet part 

After cancellation of the singularities thanks to the "+" operation the kernels are ready for 
practical usage. At a; < —^ the non-singlet part of the evolution equation reads: 



dH^s{x,^,Q^) 2as{Q' 



dlnQ^ 



37T 



dy 



"--1 



2e -x^-y^ 
{y - x){y^ - ^^) 



(HNs{y,i,Q^)-HNs{x,i,Q' 



x-i 



+ i7^s(x,e,Q')k + 21n(l + x) + — ^ln((e-a;)(l-0) (27) 



x + i 



2i 

ln((-e-a:)(l+0) 



At -^ < a; < ^: 

dHr,s{x,^,Q') a^iQ' 



cilng2 



371 



dy^77—vT7-. — -:^{HNs{y,^,Q ) -Hns{x,^,Q , 



+ / dy 



i{y^i){y-x) 

[x-i){y-x-2i) 

^{y-Oiy-x) 



H^s{y,^,Q')-H^s{x,^,Q')) (28 



+ H^six, e, Q') (3 + 2 In ^^ + ^ ln(e ' ^) ' ^ Hx + 



x>^: 



dHNs{x,^,Q^) 2a,(g2 



dlnQ^ 



3n 



dy 



x^ + y^- 2e 

(y-x)(y2-^2) 



{H^s{y,^,Q')-HNs{x,^,Q' 



X — ^ 



+ H^s{x,^,Q'){- + 2\n{l-x) + —:^\ni{x- 0(1 + 0) (29) 



x + ^ 
2e 



2e 

In ((x + 0(1-0) 



A. 2 Singlet quark-quark part 

In the DGLAP region singlet quark-quark evolution is governed by the same kernel as the 
non-singlet (Eq ^^. At < x < ^ 



dQHQ{x,^,Q^) a,{Q' 



cilnQ2 



Stt 



2x }HQ{y,^,Q')^ ,^ }HQ{y,^,Q')-HQ{x,^,Q') 



e 



^ + y 



-dy + Ax 



y2 _ ^2 



dy 



+2- 



e-x jHQ{y,^,Q' 



e { e-y' 



I , ,.,2 2^ f HQ{y,^,Q')-HQ{x,^,Q'] 
-ydy -A{i -x) — :^-— 

I (r-r)(r-a;^) 

+ 2/7Q(x,e,Q')(ln[^ + 21n| + ^^ 



ydy (30) 



A. 3 Quark-gluon part 

In this part and also in the gluon-gluon part of the kernel a pole at y = ^ must be regularized. 
For this, subtraction of the Hc{^, ^, Q^) is done. If Hg{x, C,, Q"^) is subtracted instead, this lead 
to larger of inaccuracy of the Simpson's method since the function 



{HG{y.i.Q^)-HG{x,i,Q^))/{y-i) 



varies very rapidly around y = x. 
At < a; < ^: 



dGHQ{x,^,Q') 
rflnQ2 



^f- 



2ne 



e + ex - 2ex' - {e - ix)y 



{e-y 



2A2 






{HG{y,^,Q')-HG{^,^,Q'))dy 

e — x'^ i + x^ 



X 



+ 



In 



i + x 2i i-x^ 



x>C- 

dGHQ{x,^,Q^] 
dlnQ^ 



Uf- 



«s(g' 



r 1 



27r 



1/2 - 2xy + 2a;2 - ^ 



{e~y 



2\2 



HG{y,i.Q')-HG{i.i.Q'))dy 



-HGi,,,,Q^)^'''^'-^'^'^^-'^ 



2e 



In 






(e + a;)(l-0 e(e-l). 



A. 4 Gluon-quark part 

< X < ,^: 

dQHGix,^,Q') 2as{Q') \ {^ - xf J HQ{y,tQ') - HQit^Q') 



rflnQ2 



Stt 



e 



2 ,,2 



e-y 



ydy 



+ 



e + 2iy - X- 

ei^v) 



'-HQ{y,i,Q^)dy + HQ{i,i,Q^)^^ ""^ 



In- 



e 



2e e - X 



x>0 



dQHGix, e, Q') _ 2asiQ^) l-e-x^- 2l/' + 2xy 



d\nQ^ 



37r 



e-y^ 



-HQiy,^,Q')dy. 



(31) 



(32) 



(33) 



(34) 



A. 5 Gluon-gluon part 

For the pure gluonic part it was found that simple subtraction of the GPD value ai y = x 
does not provide high enough precision. To overcome the difficulty, a partial fractioning was 
introduced so that integrand in each of the integrals had only one pole {at y = x or y = C,)- 
This brought formal complication into the evolution kernel. However, practical numerical 
calculations became more stable and precise. 
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< X < ^: 



dGHG{x,^,Q^) 3a,(g2 



rflnQ2 



4:71 



2 HG{y,^,Q )dy 



i{i + yy 



+ 






+ 4 / ^'^.+ y/'^!^ (Haiy, ^, Q') - Ha{x, ^, Q'\ 



HG{y,^,Q')-HG{^,^,Q' 



e { ie-y'r 



- 4- 



2 I ^„,2 



^ — X f S, X + '2.iy + xy 



e + xy (e-y 



2^2 



\i + x) 



2 I ^,,2 



1 /• e'x + 2^1/' + xi/ 



e-y' 



{HG{y,^,Q')-HG{^,^,Q') 
{HG{y,^,Q')-HG{^,^,Q') 



^j^^xfJ ~ — ~ [Hciy.i.Q )-Hg{x,^,Q 



y2 _ r£l 



dy 
dy 
dy 
dy 
dy 



- ^TT^(HGi^,^.Q')-HGix,^,Q': 

^ + x ^ 



X 



4.-^ (Hg{^, e, Q') - Hg{x, e, Q')) U{^ + x) In ^ + "^ 



'i^ + xy 



^-x 2 4-x^ 
a;(2^ + x) 



^-x 



HG{^,^,Q''){2{^-x)^ + l-{2^ + x){^-xr\n 



2\ I nlc ^N J^ _|_ _!_('ot I ^\/t ^n2 i._ ^^ X 



+ 2HG{x,^,Q^){\n 



l-x' 



i-x 

11 nj\ 



[l + if i + l 6 9 



x>{: 



(35) 



dGHG{x,^,Q^) SasiQ' 



dlnQ^ 



TT 



y^ + x^ + xy — X^ 

i^-ym+yy 



{HGiy,^,Q')-HGi^,^,Q'))dy 



^ r HGiy,^,Q')-HGix,^,Q') ^^ 
J y — X 

9 9 






(36) 



2e^x-x3 ^^ (e-a:)(e + l) 
4^ ""i^ + xy^-l] 



2e 



{Hg{^,^,Q')-Hg{x,^,Q 
- HG{x,i,Q 



+ 



X 



{i+x){i + i) e-i i+x 



^.{^x'+e){i-x) , / 1 



2e(e - 1) 



2^^ # '^-^'^'"^ 



1 1 - f 2 /I 

--lllT^ ^ + 



X 



2 {l-xf \2i'^ 4^ 



" '^•^ ' ^^ ^ a; - ^ 12 18 
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A. 6 Spin-dependent singlet quark-quark part 

< X <^: 



dQHQ{x,^,Q'') 2a,(Q2 



dlnQ2 



Stt 



2 I „,2 



2^y + x' + y 



{HQ{y,^,Q')-HQ{x,^,Q'))dy 



(e - ^) / ..?^''.twL^''L K(?/. e. Q') - ^q(^' e- Q')) dy (37) 



(^2_^2)(^2_3,2) 



+ ffQ(x,e,Q') In 



l-x' 



+ 



C, — X C, + X 3 



In 



+ 



(e + x)(e + i) 2e e-x 2 



For X > ^ the evolution kernel is given by the same expression as the non-singlet spin- 
independent one (Eq. 12 



A. 7 Spin-dependent quark-gluon part 

< X < ^: 

dGHQ{x,^,Q') 



dlng2 



nf- 



27r 



X 



e{^+x) 



HG{y,^,Q' 

i^+yy 



■dy 



X > ^■. 



dGHQ{x,i,Q^) as{Q' 

- = -%- 



d\nQ^ 



27r 



i +y — 2xy 

(e - y'y 



{HGiy,^,Q')-HGi^,^,Q'))dy 

1 + x 



+ 



e-i 



HG^Li.Q' 



A. 8 Spin-dependent gluon-quark part 

Q <x <i: 



dQHG{x,^,Q^) 2a,(Q2 



dlnQ^ 



3n 



2./^^|^..-ii^/^^i^M!),,, 



e { e-y' 



X > ^: 



dQHGix, e, Q') _ 2a,(g2) . ^2 + x2 - 2xy 



d\nQ^ 



Svr 



e-y^ 



^Q(l/,e,Q')^l/. 



(3^ 



(39) 



(40) 



(41) 



A. 9 Spin-dependent gluon-gluon part 

As it was in the case of the spin-independent gluon-gluon kernel, the spin- dependent gluon- 
gluon kernel requires partial fractioning of the denominator so that each of the integrals would 
have only only one pole, which can be easily handled by subtraction. 
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< X < ^: 



cilnQ2 



271 



2- I ^]^^%lT,.f {Hc{y,i.Q')-Ho{..(.Q-))dy 

,2 c2 



^-x ry^-i^-2x'- 2^a; 



i + x 



{e-y 



2^2 



HG{y.i.Q')~HG{i.i,Q^))ydy 



x>C- 



+ 



+ 



i^ + xf 



i^ + x) 



y^-i^- 2x' - 2ix 



HG{y.i.Q')~HG{i.i.Q'))ydy (42) 



y^ -i^ - 2x^ - 2ix 



y2 _ ^2 



HG{y,^,Q')-HG{x,^,Q^))ydy 



- {Hg{^, e, Q') - Hg{x, e, Q')) i Jf. . + "^ 



2 t2 

+ ln- ^ 



Ve(e+a:) {i+xf e-x' 

~. . ^2^ f 2x(l-a;) 2x^ , x2(l-a;)V 



ciGi^G(a:,e,Q') 3a,(g2 



rflng2 



27r 



2 ^ HGiy,^,Q')-HGix,^,Q') ^^^^ 



y-x 



+2 
+ HG{x,i,Q 



x{e + y') - i/(a:' + y' - e) 



(e - y 



2^2 



HG{y,i,Q^)-HG{i,i.Q^))dy (43) 



2, ,2a; -2 , (1 - x)^ 



]^ _ £2 ' "" 1 — £2 / ' ~ ^^'-'V'^' -5 5 -ts j I £2 _ ]^ 



+ 2ifG(e,e,Q' 



+ (/7g(|/, e, Q') - ^G(e, e, Q')) (27^ - 2 + 



In 
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